OpenEO Platform - User Story 3

On-demand processing of Sentinel-1 data

• CARD4L compliant data
• With custom parametrization
• Sigma to Gamma ratio
• Timeseries comparison with different pre-processing toolchains

from eo_utils import *

Step 1 - Area of Interest selection

Run the next cell, click on the rectangle button and then draw the area over the map like in the following example image

center = [46.49, 11.35]
zoom = 12

eoMap = openeoMap(center,zoom)
eoMap.map

Get the bounding box from the previous map.

bbox = eoMap.getBbox()
print('west',bbox[0],'\neast',bbox[2],'\nsouth',bbox[1],'\nnorth',bbox[3])

west 11.279467 
east 11.397167 
south 46.466356 
north 46.525914

Step 1 - openEO graph creation for CARD4L ARD Sentinel-1 processed with SentinelHub

1. Open the connection with openEO back-end. Please use the openeo-auth tool to configure credentials. https://open-eo.github.io/openeo-python-client/auth.html#config-files-and-openeo-auth-helper-tool

vito_backend_url = "https://openeo.vito.be"
vito_conn = openeo.connect(vito_backend_url).authenticate_basic()

1. Define range of interest in time for Sentinel-1

collection      = 'SENTINEL1_GRD'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2021-01-06", "2021-01-08"]
bands           = ["VV","VH"]

s1 = vito_conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)

1. Apply the CARD4L compliant SAR processing.

s1bs = s1.ard_normalized_radar_backscatter(elevation_model="COPERNICUS_30")

1. Save process as final step, here we use NetCDF as output format.

s1bs_netcdf = s1bs.save_result(format="NetCDF")

1. Now we will create a batch job and start it. More info about batch jobs can be found here:

   https://openeo.org/documentation/1.0/glossary.html#data-processing-modes

   https://open-eo.github.io/openeo-python-client/batch_jobs.html

job_vito = vito_conn.create_job(s1bs_netcdf.graph)
job_id_vito = job_vito.job_id
if job_id_vito:
    print("Batch job created with id: ",job_id_vito)
    job_vito.start_job()
else:
    print("Error! Job ID is None")

Batch job created with id:  48fecbb1-a3ec-4031-b8d2-6e1a034459b1

If our area of interest is small, we can also do a direct request, but this will not return the CARD4L json metadata. Note that this step automatically adds the save_result process at the end based on the output format we choose.

%time s1bs.download("./data/US3/card4l_out.nc",format="NetCDF")

CPU times: user 820 ms, sys: 444 ms, total: 1.26 s
Wall time: 1min 30s

We can get a description of the job and check its status.

job_vito = vito_conn.job(job_id_vito)
job_description = job_vito.describe_job()
print("Batch job with id: ",job_id_vito, ' is ',job_description['status'])

Batch job with id:  48fecbb1-a3ec-4031-b8d2-6e1a034459b1  is  running

results = job_vito.get_results()
results

results.get_metadata()

{'assets': {'out': {'card4l:nodata': 0.0,
   'eo:bands': [{'center_wavelength': None, 'name': 'VV'},
    {'center_wavelength': None, 'name': 'VH'},
    {'center_wavelength': None, 'name': 'mask'},
    {'center_wavelength': None, 'name': 'local_incidence_angle'}],
   'href': 'https://openeo.vito.be/openeo/1.0/jobs/6d266445-6a10-46b3-a986-539365cc999f/results/out',
   'type': 'application/x-netcdf'},
  's1_rtc_02FB1E_N46E011_2021_01_07_MULTIBAND.tif': {'href': 'https://openeo.vito.be/openeo/1.0/jobs/6d266445-6a10-46b3-a986-539365cc999f/results/s1_rtc_02FB1E_N46E011_2021_01_07_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_02FB1E_N46E011_2021_01_07_metadata.json': {'href': 'https://openeo.vito.be/openeo/1.0/jobs/6d266445-6a10-46b3-a986-539365cc999f/results/s1_rtc_02FB1E_N46E011_2021_01_07_metadata.json',
   'type': 'application/json'},
  's1_rtc_04393F_N46E011_2021_01_08_MULTIBAND.tif': {'href': 'https://openeo.vito.be/openeo/1.0/jobs/6d266445-6a10-46b3-a986-539365cc999f/results/s1_rtc_04393F_N46E011_2021_01_08_MULTIBAND.tif',
   'type': 'image/tiff; application=geotiff'},
  's1_rtc_04393F_N46E011_2021_01_08_metadata.json': {'href': 'https://openeo.vito.be/openeo/1.0/jobs/6d266445-6a10-46b3-a986-539365cc999f/results/s1_rtc_04393F_N46E011_2021_01_08_metadata.json',
   'type': 'application/json'}},
 'bbox': [11.255493, 46.465358, 11.377716, 46.528698],
 'geometry': {'coordinates': [[[11.255493, 46.465358],
    [11.255493, 46.528698],
    [11.377716, 46.528698],
    [11.377716, 46.465358],
    [11.255493, 46.465358]]],
  'type': 'Polygon'},
 'id': '6d266445-6a10-46b3-a986-539365cc999f',
 'links': [{'href': 'https://openeo.vito.be/openeo/1.0/jobs/6d266445-6a10-46b3-a986-539365cc999f/results',
   'rel': 'self',
   'type': 'application/json'},
  {'href': 'http://ceos.org/ard/files/PFS/SR/v5.0/CARD4L_Product_Family_Specification_Surface_Reflectance-v5.0.pdf',
   'rel': 'card4l-document',
   'type': 'application/pdf'}],
 'properties': {'created': '2021-03-19T09:14:45Z',
  'datetime': None,
  'end_datetime': '2021-01-08T00:00:00Z',
  'instruments': ['c-sar'],
  'processing:lineage': {'process_graph': {'ardnormalizedradarbackscatter1': {'arguments': {'contributing_area': False,
      'data': {'from_node': 'loadcollection1'},
      'elevation_model': 'COPERNICUS_30',
      'ellipsoid_incidence_angle': False,
      'noise_removal': True},
     'process_id': 'ard_normalized_radar_backscatter'},
    'loadcollection1': {'arguments': {'bands': ['VV', 'VH'],
      'id': 'SENTINEL1_GRD',
      'spatial_extent': {'crs': 'EPSG:4326',
       'east': 11.377716,
       'north': 46.528698,
       'south': 46.465358,
       'west': 11.255493},
      'temporal_extent': ['2021-01-06', '2021-01-08']},
     'process_id': 'load_collection'},
    'saveresult1': {'arguments': {'data': {'from_node': 'ardnormalizedradarbackscatter1'},
      'format': 'NetCDF',
      'options': {}},
     'process_id': 'save_result',
     'result': True}}},
  'proj:epsg': 3785,
  'start_datetime': '2021-01-06T00:00:00Z',
  'updated': '2021-03-19T09:23:01Z'},
 'stac_extensions': ['processing', 'card4l-eo', 'eo', 'projection'],
 'stac_version': '0.9.0',
 'type': 'Feature'}

Download complete result with metadata and source files

CARD4L results contain STAC metadata and the requested image.

We can simply download everything, for inspection. Please note: this will download also the original S1 GRD files, which can take a lot of time for a big area or a timeseries.

# results.download_files()

The downloaded data can be opened, but these are fairly large files, making visualization more difficult.

The foreseen way of interacting with it is to use openEO to further process the dataset into a more manageable result.

# sar = xr.open_rasterio('s1_rtc_02F8D2_N46E011_2021_01_02_MULTIBAND.tif')
# sar

<xarray.DataArray (band: 4, y: 5000, x: 5000)>
[100000000 values with dtype=float32]
Coordinates:
  * band     (band) int64 1 2 3 4
  * y        (y) float64 47.0 47.0 47.0 47.0 47.0 ... 46.0 46.0 46.0 46.0 46.0
  * x        (x) float64 11.0 11.0 11.0 11.0 11.0 ... 12.0 12.0 12.0 12.0 12.0
Attributes:
    transform:               (0.0002, 0.0, 11.0, 0.0, -0.0002, 47.0)
    crs:                     +init=epsg:4326
    res:                     (0.0002, 0.0002)
    is_tiled:                1
    nodatavals:              (nan, nan, nan, nan)
    scales:                  (1.0, 1.0, 1.0, 1.0)
    offsets:                 (0.0, 0.0, 0.0, 0.0)
    AREA_OR_POINT:           Area
    TIFFTAG_RESOLUTIONUNIT:  1 (unitless)
    TIFFTAG_XRESOLUTION:     1
    TIFFTAG_YRESOLUTION:     1

Output visualization

S1_ard = xr.open_dataset("./data/US3/card4l_out.nc")
S1_ard

<xarray.Dataset>
Dimensions:                (t: 2, x: 843, y: 843)
Coordinates:
  * t                      (t) datetime64[ns] 2021-01-07 2021-01-08
  * x                      (x) float64 6.757e+05 6.757e+05 ... 6.842e+05
  * y                      (y) float64 5.147e+06 5.147e+06 ... 5.156e+06
Data variables:
    VH                     (t, y, x) float32 ...
    VV                     (t, y, x) float32 ...
    local_incidence_angle  (t, y, x) float32 ...
    mask                   (t, y, x) float32 ...
Attributes:
    nodata:   0.0
    crs:      +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 

fig, (ax1, ax2) = plt.subplots(1,2,figsize=(20,20))
ax1.imshow(np.flipud(10* np.log10(S1_ard.VV[1].values)),cmap='Greys_r',vmin=-30,vmax=30)
ax1.set_title('VV gamma0')
ax2.imshow(np.flipud(10* np.log10(S1_ard.VH[1].values)),cmap='Greys_r',vmin=-30,vmax=30)
ax2.set_title('VH gamma0')

Text(0.5, 1.0, 'VH gamma0')

If we don't need CARDL compliant data, we can also use the sar_backscatter process:

collection      = 'SENTINEL1_GRD'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2021-01-01", "2021-01-05"]
bands           = ["VV","VH"]

s1 = vito_conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)

Sigma0:

s1_sigma = s1.sar_backscatter(coefficient='sigma0-ellipsoid').apply(lambda x: 10 * x.log(base=10))

%time s1_sigma.download("./data/US3/s1_sigma.nc",format="NetCDF")

CPU times: user 180 ms, sys: 224 ms, total: 404 ms
Wall time: 30 s

Gamma0:

s1_gamma = s1.sar_backscatter(coefficient='gamma0-ellipsoid').apply(lambda x: 10 * x.log(base=10))

%time s1_gamma.download("./data/US3/s1_gamma.nc",format="NetCDF")

CPU times: user 76 ms, sys: 112 ms, total: 188 ms
Wall time: 31.3 s

s1_sigma0_ = xr.open_dataset('./data/US3/s1_sigma.nc')
s1_gamma0_ = xr.open_dataset('./data/US3/s1_gamma.nc')

fig, (ax1, ax2) = plt.subplots(1,2,figsize=(20,20))

ax1.imshow(np.flipud(s1_sigma0_.VV[1].values),cmap='Greys_r',vmin=-30,vmax=30)
ax1.set_title('Sigma VV')
ax2.imshow(np.flipud(s1_sigma0_.VH[1].values),cmap='Greys_r',vmin=-30,vmax=30)
ax2.set_title('Sigma VH')

Text(0.5, 1.0, 'Sigma VH')

fig, (ax1, ax2) = plt.subplots(1,2,figsize=(20,20))

ax1.imshow(np.flipud(s1_gamma0_.VV[1].values),cmap='Greys_r',vmin=-30,vmax=30)
ax1.set_title('Gamma VV')
ax2.imshow(np.flipud(s1_gamma0_.VH[1].values),cmap='Greys_r',vmin=-30,vmax=30)
ax2.set_title('Gamma VH')

Text(0.5, 1.0, 'Gamma VH')

Sigma To Gamma Ratio

fig, (ax1, ax2) = plt.subplots(1,2,figsize=(20,20))

ax1.imshow(np.flipud(s1_sigma0_.VV[1].values/s1_gamma0_.VV[1].values),cmap='Greys_r',vmin=0,vmax=1)
ax1.set_title('Sigma/Gamma VV')
ax2.imshow(np.flipud(s1_sigma0_.VH[1].values/s1_gamma0_.VH[1].values),cmap='Greys_r',vmin=0,vmax=1)
ax2.set_title('Sigma/Gamma VH')

Text(0.5, 1.0, 'Sigma/Gamma VH')

Timeseries comparison

Next to requesting CARD4L backscatter data on SentinelHub, we can also generate Sigma0 backscatter on CreoDIAS.

backend_url = "https://openeo.creo.vito.be"
creo_conn = openeo.connect(backend_url).authenticate_basic()

W, S = 510000, 5680000
bbox = {
    "west": W, "east": W + 32 * 10,
    "south": S, "north": S + 32 * 10,
    "crs": 32631
}

dates = ("2020-05-06T00:00:00", "2020-09-30T00:00:00")

def backscatter(connection):
    return (connection.load_collection("SENTINEL1_GRD")
    .filter_bbox(**bbox)
    .filter_temporal(dates)
    .filter_bands(["VH", "VV"])
    .sar_backscatter(coefficient="sigma0-ellipsoid")
    .apply(lambda x: 10 * x.log(base=10)))

%%time
backscatter(creo_conn).download("./data/US3/creo-series.nc", format="NetCDF")

CPU times: user 48 ms, sys: 8 ms, total: 56 ms
Wall time: 2min 27s

%%time
backscatter(vito_conn).download("./data/US3/shub-series.nc", format="NetCDF")

CPU times: user 36 ms, sys: 12 ms, total: 48 ms
Wall time: 20.4 s

%%time
asc = (vito_conn.load_collection("S1_GRD_SIGMA0_ASCENDING").filter_bbox(**bbox)
    .filter_temporal(dates)
    .filter_bands(["VH", "VV"]))
desc = (vito_conn.load_collection("S1_GRD_SIGMA0_DESCENDING").filter_bbox(**bbox)
    .filter_temporal(dates)
    .filter_bands(["VH", "VV"]))
desc.merge_cubes(asc,overlap_resolver="max").apply(lambda x: 10 * x.log(base=10)).download("./data/US3/snap-series.nc", format="NetCDF") 

CPU times: user 20 ms, sys: 8 ms, total: 28 ms
Wall time: 12 s

creo_ts = xr.open_dataset("./data/US3/creo-series.nc",engine="h5netcdf")
creo_ts

<xarray.Dataset>
Dimensions:  (t: 98, x: 32, y: 32)
Coordinates:
  * t        (t) datetime64[ns] 2020-05-07 2020-05-08 ... 2020-09-28 2020-09-29
  * x        (x) float64 5.1e+05 5.1e+05 5.1e+05 ... 5.103e+05 5.103e+05
  * y        (y) float64 5.68e+06 5.68e+06 5.68e+06 ... 5.68e+06 5.68e+06
Data variables:
    VH       (t, y, x) float32 ...
    VV       (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs 
    nodata:   nan

shub_ts = xr.open_dataset("./data/US3/shub-series.nc",engine="h5netcdf")
shub_ts

<xarray.Dataset>
Dimensions:  (t: 98, x: 32, y: 32)
Coordinates:
  * t        (t) datetime64[ns] 2020-05-07 2020-05-08 ... 2020-09-28 2020-09-29
  * x        (x) float64 5.1e+05 5.1e+05 5.1e+05 ... 5.103e+05 5.103e+05
  * y        (y) float64 5.68e+06 5.68e+06 5.68e+06 ... 5.68e+06 5.68e+06
Data variables:
    VH       (t, y, x) float32 ...
    VV       (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs 
    nodata:   nan

snap_ts = xr.open_dataset("./data/US3/snap-series.nc",engine="h5netcdf")
snap_ts

<xarray.Dataset>
Dimensions:  (t: 98, x: 32, y: 32)
Coordinates:
  * t        (t) datetime64[ns] 2020-05-07 2020-05-08 ... 2020-09-28 2020-09-29
  * x        (x) float64 5.1e+05 5.1e+05 5.1e+05 ... 5.103e+05 5.103e+05
  * y        (y) float64 5.68e+06 5.68e+06 5.68e+06 ... 5.68e+06 5.68e+06
Data variables:
    VH       (t, y, x) float32 ...
    VV       (t, y, x) float32 ...
Attributes:
    crs:      +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs 
    nodata:   nan

Compute the mean value, for each timestamp. OpenEO can also do this for you.

xr.merge([creo_ts.rename({"VV": "VV_ORFEO", "VH":"VH_ORFEO"}), shub_ts.rename({"VV": "VV_SHUB", "VH":"VH_SHUB"}), snap_ts.rename({"VV": "VV_SNAP", "VH":"VH_SNAP"})]).mean(dim=['x','y']).hvplot(width=1000,height=600)

Step 2 - openEO graph creation for CARD4L ARD Sentinel-1 processed with SNAP

1. Open the connection with openEO back-end. Please use the openeo-auth tool to configure credentials. https://open-eo.github.io/openeo-python-client/auth.html#config-files-and-openeo-auth-helper-tool

eodc_backend_url = "https://openeo.eodc.eu/v1.0"
eodc_conn = openeo.connect(eodc_backend_url).authenticate_basic()

Explore the pre-computed data and find the available dates over our area of interest:

itemIter = eodc_conn.collection_items(
    "gamma0_sentinel_1_dv",
    temporal_extent = ["2017-03-14", "2017-03-21"],
    spatial_extent = bbox
)

next(itemIter)

1. Define range of interest in time for Sentinel-1

collection      = 's1a_csar_grdh_iw'
spatial_extent  = {'west':bbox[0],'east':bbox[2],'south':bbox[1],'north':bbox[3],'crs':'EPSG:4326'}
temporal_extent = ["2017-03-19", "2017-03-20"]
bands           = ["VV+VH"]

s1 = eodc_conn.load_collection(collection,spatial_extent=spatial_extent,bands=bands,temporal_extent=temporal_extent)

1. Apply the CARDL compliant SAR processing.

# s1bs = s1.ard_normalized_radar_backscatter(elevation_model="COPERNICUS_30")
s1bs = s1.sar_backscatter()

1. Save process as final step, here we use GeoTIFF as output format.

s1bs_tif = s1bs.save_result(format="GTiff")

1. Now we will create a batch job and start it. More info about batch jobs can be found here:

   https://openeo.org/documentation/1.0/glossary.html#data-processing-modes

   https://open-eo.github.io/openeo-python-client/batch_jobs.html

job_eodc = eodc_conn.create_job(s1bs_tif.graph)
job_id_eodc = job_eodc.job_id
if job_id_eodc:
    print("Batch job created with id: ",job_id_eodc)
    job_eodc.start_job()
else:
    print("Error! Job ID is None")

Batch job created with id:  jb-b1fa0f26-0cd2-4b8d-8c3b-7db4a6cc0726

job_eodc = eodc_conn.job(job_id_eodc)
job_description = job_eodc.describe_job()
print("Batch job with id: ",job_id_eodc, ' is ',job_description['status'])

Batch job with id:  jb-b1fa0f26-0cd2-4b8d-8c3b-7db4a6cc0726  is  running

results = job_eodc.get_results()
results

results.download_files()

[PosixPath('/home/jovyan/SRR1_updated/SRR1_notebooks/S1A__IW___A_20170319T170640_DVH_NR_Orb_TF_TC_EU010M_E048N014T1.tif'),
 PosixPath('/home/jovyan/SRR1_updated/SRR1_notebooks/S1A__IW___A_20170319T170640_DVV_NR_Orb_TF_TC_EU010M_E048N014T1.tif')]

Output visualization

S1_ard_snap = xr.open_rasterio("S1A__IW___A_20170319T170640_DVH_NR_Orb_TF_TC_EU010M_E048N014T1.tif")

fig, ax1 = plt.subplots(1,1,figsize=(10,10))
ax1.imshow(10* np.log10(S1_ard_snap.where(S1_ard_snap>0,drop=True)[0].values),cmap='Greys_r')
ax1.set_title('VV gamma0 SNAP')

Text(0.5, 1.0, 'VV gamma0 SNAP')